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We compute the two-particle quantities relevant for superconducting correlations in the two- 
dimensional Hubbard model within the dynamical cluster approximation. In the normal state we 
identify the parameter regime in density, interaction, and second-nearest-neighbor hopping strength 
that maximizes the d„, 2 _y 2 superconducting transition temperature. We find in all cases that the 
optimal transition temperature occurs at intermediate coupling strength, and is suppressed at strong 
and weak interaction strengths. Similarly, superconducting fluctuations are strongest at intermediate 
doping and suppressed towards large doping and half-filling. We find a change in sign of the vertex 
contributions to dj;y superconductivity from repulsive near half filling to attractive at large doping, 
p-wave superconductivity is not found at the parameters we study, and s-wave contributions are 
always repulsive. For negative second-nearest-neighbor hopping the optimal transition temperature 
shifts towards the electron-doped side in opposition to the van Hove singularity which moves towards 
hole doping. We surmise that an increase of the local interaction of the electron-doped compounds 
would increase Tc- 


Understanding physical scenarios that give rise to su¬ 
perconductivity at high temperatures has been a pri¬ 
mary motivating force behind computational research of 
strongly correlated electron systems and candidate mod¬ 
els such as the 2D Hubbard model J . Only recently 
have reliable many-body methods 3j become powerful 
enough to reach temperatures low enough to cross the 
superconducting transition at intermediate interaction 
strengths [3-0 , but progress is limited by the exponential 
scaling intrinsic to all unbiased methods. Such compu¬ 
tational work has identified clearly the competition be¬ 
tween correlations that give rise to superconductivity and 
other phases such as antiferromagnetism 0, 0 and the 
pseudogap ii phenomenon within the 2D Hubbard 
model. 

Central to understanding these phases is the evalua¬ 
tion of two-particle susceptibilities and vertex functions 
at nonzero temperature, which diverge on approach to a 
continuous phase transition and may also exhibit signs of 
a transition at temperatures much larger than the transi¬ 
tion temperature, for parameters that are accessible with 
current techniques and computational power. Neverthe¬ 
less, the numerical calculation of these two-particle sus¬ 
ceptibilities requires techniques that are robust across the 
full phase diagram, can reach low temperatures, are ca¬ 
pable of providing reliable and systematically improv¬ 
able results, and are able to distinguish independent 
phases. Cluster dynamical mean field theory [9l-ll3l| pro¬ 
vides such a self consistent non-perturbative tool for sim¬ 
ulating strongly correlated electron problems. The dy¬ 
namical cluster approximation (DCA) is based on a self¬ 
energy discretization into Nc independent self-ener^ co¬ 
efficients which recover the exact limit as Nc —> oo 0113- 
m and capture much of the physics believed to be rel¬ 
evant for the superconductivity and pseudogap physics 
two-dimensional Hubbard model on clusters of size 8 and 


larger 0,0 . 

In this work, we specifically address the problem of 
optimizing the superconducting transition temperature 
in the 2D Hubbard model by analyzing wide regions of 
parameter space. We first demonstrate how the vertex 
contribution to the pairing susceptibility can be used 
as an indicator of the proximity to the superconduct¬ 
ing transition temperature, Tc- We then show that this 
quantity, as temperature is reduced, mimics the depen¬ 
dence of Tc on model parameters. This allows us to sweep 
the entirety of parameter space in density n, interaction 
strength U/t, and second-nearest neighbor hopping t'/t 
at numerically accessible T TTc, to identify regions 
of qualitatively high or low Tc, so that the maxima can 
then be targeted for a quantitative determination of the 
optimal Tc value. We mainly focus on dj, 2 _y 2 supercon¬ 
ductivity but show that dxy superconducting fluctuations 
change from attractive (at large doping) to repulsive (at 
low doping), p-wave fluctuations are always either repul¬ 
sive or zero within error bars at the system sizes, in¬ 
teraction strengths, and dopings we study, and s wave 
contributions are strongly repulsive. 

We study the single orbital Hubbard model in two di¬ 
mensions with nearest and next-nearest hopping param¬ 
eters, 

H = '^{ek- u) cl^Cka + UY^ ( 1 ) 

fc,(T i 

where p is the chemical potential, k momentum, i la¬ 
bels sites in real-space, U is the interaction, and the 
dispersion is given by Ck = —2t [cos{kx) + cos{ky)] — 
41'cos{kx) cos{ky). We operate in a formalism that al¬ 
lows for a nonzero anomalous Creen’s function in the 
superconducting state, which is defined as F(/c,r) = 
— {TrCk^{T)c-ki{0))- At T > Tc superconducting order 
will be absent but fluctuations are captured by the gener- 
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alized susceptibility, written in imaginary time in terms 
of the one- and two-particle Green’s functions as 0 (see 
Supplemental Material [T^ for definition and notations) 

Xc7'ia'2(73(T4 (^ 1^1 7 ^ 2^2 7 ^ 3^3 7 (^) 

— G2,ai. ..(74 (fcin, fc 2 T 2 , fcsTa, fc4T4) 

Gcria2 (^ 1 ^ 1 7 ^2'^2')G (^ 3^3 ; ^ 4 "^ 4 ) 
or as its Fourier transform 

Xppaa'{k,k',q)= [ f f dTidT2dT3 (3) 
Jo Jo Jo 

XXaaa'a'ikTi, {q - k')T 2 , (q - fc)r3, fc'O) 

^i{v-Uj')T2 ^-i{l'-Uj)T3 


where ui and uj' are fermionic Matsubara frequencies, 12 is 
a bosonic Matsubara frequency, cr and a' are t or spin 
labels and fc, k' and q are initial, final and transfer mo¬ 
menta respectively, and pp denotes the Fourier transform 
convention. With the difference between the aa' =tt and 
susceptibilities defined as Xppfi = XppTT “ XppTt > linear 
response theory relates XppJl to the response of a system 
to a generating superconducting field r]{k) 


dr 


SF(k', T = 0; 77 ) 


Sr]{k,T) 


7J—0 


= ® If *'.« = “) 


(4) 

where F{k', r; 77 ) is the anomalous Green’s function com¬ 
puted in the presence of an external superconducting 
field. We note that the quantity on the left-hand side 
is commonly referred to as the uniform pairing suscepti¬ 
bility l3, ifl]- 

Continuous phase transitions can be identified by the 
point in phase space where the corresponding suscepti¬ 
bility diverges. The susceptibility can then, using the 
Bethe-Salpeter equation, be separated into a ‘bare’ con¬ 
tribution 


u;u; ly 
Xopp 


{k, k', q) = -PGa{k,iu})Ga{q - k',iv - iLo')5u!u>'5kk' 

(5) 


which never diverges and a part containing an irreducible 
vertex function Fpp, 


X 


ujuj' y 

ppJl 


ik,k',q) 


= Xtpp^iK k\ q) - ^x;;ifik, k", q) 

<l)xC‘^''^ik'”, k\ q). ( 6 ) 


In order to see the origin of the divergence in Xp^^ this 
susceptibility can be expressed in matrix notation giving 


— = Xo 


(7) 


and the point of divergence of x is identified as the point 
where an eigenvalue of --^^ppjiXo crosses 1 , and the 



FIG. 1. Left panel: Superconducting critical temperature 
of the Hubbard model with nearest neighbour hopping and 
next nearest neighbour hopping t' = —O.lt for U = 6t using 
an Nc — 8 dynamical cluster approximation. Right panel: 
Pd^ 2 _ 2 8-t different temperatures with interaction strength 
U — 6t and next nearest neighbour hopping t' = —O.lt using 
an Ale = 8 cluster. 


symmetry of the eigenvector will identify the symmetry 
of the order parameter. 


In what follows we solve the Hubbard model within the 
(paramagnetic) dynamical cluster approximation which 
approximates the self-energy of the interacting model by 
a number, Nc, of ‘coarse-grained’ frequency-dependent 
but momentum-independent self-energy tiles. We pri¬ 
marily present results for an Nc = 8 cluster since this is 
the smallest DGA system that captures a clear distinc¬ 
tion between nodal and antinodal physics 2ll-l24| . Gom- 
parisons to larger and smaller Nc = A and Nc = 16 sys¬ 
tems are shown in the supplemental materials 0 . An¬ 
tiferromagnetic order is actively suppressed in our cal¬ 
culations by enforcing paramagnetic spin symmetry, and 
the presence or effect of charge order [ 2 ^ has not been 
investigated. The DGA calculation provides one- and 
two-particle cluster Green’s functions, from which we ex¬ 
tract cluster susceptibilities and, using the formalism out¬ 
lined in Ref. [i^, the coarse-grained lattice susceptibil¬ 
ities x“^(Ar, AT', Q), where K, K' and Q are cluster 
momenta. In order to analyze the angular dependence 
of the superconducting order, one typically performs a 
multipole expansion restricted to the D^h square lattice 
symmetry. |27 h 29| | Because of our limited momentum res¬ 
olution we project out and analyze the leading contribu¬ 
tion and are insensitive to higher order harmonics around 
the Fermi surface. The accessible s—, p—,dxy or d^^-y^ 
symmetries are enforced by including symmetry factors 
g{K)g{K') while summing over all initial K and final K' 
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FIG. 2. Pd 2 2 different interaction strengths as a function of carrier concentration on an eight-site cluster, for t' = 0 
(panel a), t' = —O.lt (panel b) and t' = —0.2f (panel c) at /3 = 15/t. U = At (black solid line, circles), 5t (red dotted line, 
squares), 6t (green dashed line, diamonds) and 7t (blue dash-dotted line, triangles). 


states in 


Eq. O [lilaSill, with 


9iK) 


1 s 

sin(ii:,^,) p 

sm{Ka:)sin{Ky) d^y 

cos{Kx) — cos{Ky) dx 2 -y 2 


( 8 ) 


The divergence of ^ is caused by the vertex correc¬ 
tion part ~ Xo‘^ ■ We impose a shorthand nota¬ 

tion for this quantity of interest, which we call the cor¬ 
related pairing susceptibility Pg, where g refers to the 
corresponding symmetry function defined in Eq. |S1 and 
we take this to be the summation over fermionic Mat- 
subara frequencies and momenta: 


Pg :=(X - Xo)g = X! 9 {K)g{K') 


(9) 


cjcj'KK' 


K', 0 ) - xr'°iK, K', 0)1 / 


K 


We show in the supplemental material [l^ an explicit ex¬ 
ample where the point of divergence of x coincides with 
the divergence of a single eigenvalue with d 2 ; 2 _y 2 symme¬ 
try. 

The fact that the correlated pairing susceptibility Pg 
must become large on approach to Tf. grants us additional 
insights at T > Tc, where Pg can be used as a qualitative 
measure of the proximity of the system to a transition. 
The left panel of Fig. [1] shows the critical temperature 
obtained from systematically reducing T and explicitly 
evaluating the eigenvalues of ~^^pp'fiXo to find the di¬ 
vergence of the da; 2 _y 2 susceptibility. The right panel 
contrasts this with the magnitude of Pg at much higher 
temperatures (3 = 10/t, 15/t, 20/t, and 25/t. We see Pg 
tracks Tc and shows the largest superconducting fluctua¬ 
tions approximately where Tc is highest, as also indicated 
by the vertical blue lines. The correspondence of Pg to 
Tc improves as T decreases towards Tc. 

In Fig. [2Ka) we explore Pd^ 2 _y 2 a function of par¬ 
ticle density n {n = 1 denotes half filling) in the inter¬ 
mediate interaction strengths regime C//t = 4 to 7 at 


P = 15/f « 2Tc. For the weakest interaction strength 
considered here, U = 4t, the superconducting d 2 ; 2 _y 2 
fluctuations are strongest at half filling and decrease 
rapidly towards larger hole and electron doping. At 10% 
doping, the model has been shown to be dx 2 -y 2 supercon¬ 
ducting by DCA calculations extrapolated to the thermo¬ 
dynamic limit 32| , and 8-site fluctuations have shown to 
be weaker than for the lattice model. The maximum 
of fluctuations at half filling is consistent with results 
from weak coupling theory [s^, FLEX [s^, and dia¬ 
grammatic Monte Carlo calculations in the weak coupling 
limit and is also observed in results from lattice 

quantum Monte Carlo (QMC) simulation 133 and the 
two-particle-self-consistent approximation [36|. Reduc¬ 
tion of U rapidly suppresses the strength of fluctuations. 
Pd 2 9 increases at all n as U is raised to 5t. As U is 
further raised to 6t, the strength of fluctuations increases 
away from half filling but decreases near half filling, and 
the fluctuation maximum moves to finite doping, estab¬ 
lishing a dome. The suppression at half Filing coincides 
with the establishment of a pseudogap at this interaction 
strength 2lL a nd is also seen in QMC simulation [l^ 
and TPSC 130. 13711 (though it seems to be absent in four- 
site CDMFT |3a). Simulations directly in the supercon¬ 
ducting phase p have also shown that that supercon¬ 
ductivity in this region is suppressed. Above U/t = 6.4 
the half-filled system becomes Mott insulating [22j and 
superconducting fluctuations are further suppressed (but 
remain nonzero), while their maximum strength moves 
to higher doping, giving the appearance of a dome struc¬ 
ture centered at doping, 5 ~ ±1/8 for U/t ^ 8. As the 
interaction strength is further increased, fluctuations are 
suppressed and quickly decay, in qualitative agreement 
with simulations of the t — J model [3^ and Hubbard 
NLCE calculations (2^ . 

Figure [3] expands further upon the data of Fig. [2l in¬ 
cluding additional data points at intermediate interaction 
values, as a false color contour plot at = 0 in Fig.[3Ka). 
The plot clearly shows the intermediate interaction re¬ 
gion most conducive to superconductivity. The point of 
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maximum susceptibility which occurs at jg 

marked as + and occurs at U/t = 6, n = 0.92 for the 
eight-site cluster. A wide area in the vicinity of this 
point exhibits fluctuation within 10% of the maximum 
value, showing that d^ 2 _y 2 superconducting fluctuation 
is a robust feature of the model. Finite size effects change 
the precise location and general strength of the fluctua¬ 
tions (see Supplemental Material 0 ) but not the overall 
shape. Long-range antiferromagnetism may preempt the 
superconducting phase near half filling; see e.g. Ref. 0. 

Next-nearest neighbor hopping, shown in Figs. HKb), 
He) and Fig. |31 has a profound effect on fluctua¬ 

tions. As the interaction strength is raised, a pronounced 
particle hole asymmetry appears for t'/t = —0.1 (panel 
(b)) that increases superconducting fluctuations on the 
electron doped side (n > 1) while suppressing them on 
the hole doped side. Increasing t' to —0.2t (Fig. H^)) 
leads to a further enhancement of fluctuations on the 
electron doped side and increased suppression on the hole 
doped side near half filling. This behavior seems to be 
unrelated to any feature in the single particle density 
of states which has a van Hove maximum on the hole- 
doped side. Rather, we attribute it to the establishment 
of a pseudogap on the hole doped side, which is absent 
on the electron doped side j22| . and which is known to 
rapidly suppress critical temperature near half filling Q . 
The magnitude of fluctuation at the electron doped side 
(and outside of the pseudogap region at the hole-doped 
side) is not significantly changed, suggesting (in agree¬ 
ment with ED and DMRG simulations on t — J ladders 
[s^ li^ and NCA results on 2 x 2 clusters) that the f 
trends observed in real materials are not captured by the 
single band Hubbard model [4l|. We find that further 
increase of t' continues this trend and reduces the overall 
susceptibility to d^^-y^ superconductivity. 

Our results suggest that the low-energy effective mod¬ 
els of high Tc compounds do not just differ by t', but also 
by their on-site interactions U. As the electron-doped 
compounds have a much lower critical temperature than 
the hole doped ones, we surmise that they are not local¬ 
ized at the point in phase space that yields the highest 
Tc, and that an increase of U would rapidly increase the 
critical temperature. 

Finally, we establish the absence of high-temperature 
superconducting fluctuations in other symmetry channels 
by considering g{k)g{k') factors with alternate symmetry 
in the two-particle representation of the susceptibility. 
We plot results for t'/t = 0, -0.1, and -0.2 in Figs.HKS'^^c) 
at U/t = 6, for dxy and p-wave symmetry and include 
dx' 2 -y 2 for reference (also shown in Fig. [2|). 

In the large doping weak coupling regime, dx 2 -y 2 su¬ 
perconductivity is preempted by dxy superconductivity 
4^ 4^ . This is also found in RPA calculations 3 
and diagrammatic QMC calculations jsij. In contrast, 
the vertex contribution to dxy superconductivity is re¬ 
pulsive near half filling, consistent with early QMC cal¬ 



n 


FIG. 3. Contour plots for Pd 2 _ 2 in space of interaction 
strength and carrier concentration on an 8-site cluster at 
/3 = 15/t. Top panel: t' = 0. Middle panel: t' = —O.lt. 
Bottom panel: t' = -0.2t. occurs at ) = 

(5.5,0.95 and 1.05), (6,1.03), (6,1.01) respectively, marked 
by a -I- symbol. 
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FIG. 4. The correlated pairing susceptibility Pg in different 
symmetry channels with interaction strength U = 6t, at jdt = 
15, using 8-site cluster. Panel(a): t' = 0; panel(b): t' = —O.lt; 
panel(c): t' = —0.2t 


culations 0- Figure 0] shows how it changes sign for 
larger doping and eventually becomes the dominant con¬ 
tribution. 

As U is raised in the dilute (n —)■ 0) limit, dxy or¬ 
der is replaced immediately by p-wave superconductivity 
^ , [ 34 1. Third order perturbative calculations also 


find a large range of p-wave stability (but no dxy) in the 
large doping regime at 1/ = 6, and DCA calculations sim¬ 
ilarly found dominant p-wave contributions [0. Within 
our calculations, p-wave contributions to the vertex are 
zero within errors in the entire range of phase space, ex- 
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cept near half filling, where they are repulsive. Our data 
are consistent with Ref. ^ on the level of the suscep¬ 
tibility, but we find that the dominant contribution ob¬ 
served in that work is carried by xoj not the vertex part. 
Whether a DCA simulation could find dominant p-wave 
contributions to the vertex at smaller U, lower T, or on 
larger systems is an open question. The highest critical 
temperature of any non-(ia, 2 _y 2 superconductivity is far 
below the T examined in this work. 

Over the entire phase space, s-wave superconductivity 
(not plotted in Fig. |4|) is strongly repulsive, consistent 
with QMC calculations [l3, 47, 48|. At t' = —0.2 and in 
the dilute limit, weak coupling and RPA results suggest 


a favored dxy symmetry [4^ [43, l49j , consistent with our 
results at larger U/t and high temperatures. 

In summary, we have identified the regions in param¬ 
eter space that give optimal superconducting transition 
temperatures, using a formalism based on two-particle 
simulations at temperatures much higher than Tc, We 
have explored the susceptibility of the Hubbard model 
towards superconducting order over the entirety of the 
phase diagram. 

We find that both weak and strong interaction regimes, 
as well as low doping and half filled regimes, are nonop- 
timal for superconducting fluctuations, but that there is 
a large region that is very conducive to superconductiv¬ 
ity. For t' < 0 we find a shift of the optimal supercon¬ 
ducting features to the electron-doped side of the phase 
diagram, due to the establishment of a competing pseu¬ 
dogap on the hole-doped side. As actual electron-doped 
compounds have a lower Tc than the hole-doped ones, we 
surmise that a rapid increase of Tc could be achieved by 
changing the effective on-site interaction. 

By examining alternate order symmetries and t' < 0 
we show susceptibility towards d^y but not p-wave super¬ 
conductivity in the strongly hole doped n ^ 0 (dilute) 
limit. We emphasize that transitions to those symme¬ 
tries happen at temperatures much lower than the T we 
have examined here. 
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